Understanding fragility in supercooled Lennard-Jones mixtures. 
II. Potential energy surface 
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We numerically investigated the connection between isobaric fragility and the properties of high- 
order stationary points of the potential energy surface in different supercooled Lennard-Jones mix- 
tures. The increase of effective activation energies upon supercooling appears to be driven by the 
increase of average potential energy barriers measured by the energy dependence of the fraction of 
unstable modes. Such an increase is sharper, the more fragile is the mixture. Correlations between 
fragility and other properties of high-order stationary points, including the vibrational density of 
states and the localization features of unstable modes, are also discussed. 

PACS numbers: 61.43.Fs, 61.20.Lc, 64.70.Pf, 61.20.Ja 



I. INTRODUCTION 

The study of the Potential Energy Surface (PES), or 
energy landscape, of supercooled liquids and glasses is 
of fundamental importance for understanding thermody- 
namical and dynamical properties in these systems 1^*2, 
Since the pioneering work of Stillinger and Weber f2 a 
growing body of data, coming from numerical simula- 
tions, has provided a detailed description of the en- 
ergy landscape explored by supercooled liquids. As 
the temperature is lowered toward the glass transition, 
progressively deeper regions of the PES are visited)^ 
where the basins of attraction of groups of local minima 
(metabasins)^ act as traps for the system in configuration 
space^i^ii^ii^ and slow down the liquid dynamics. Accord- 
ingly, viscosity and structural relaxation times show a 
dramatic increase upon supercooling. In the so called 
fragile glass-formers, such an increase is faster than Ar- 
rhenius (or super- Arrhenius) , as opposed to the behavior 
of strong glass-formers, in which the temperature depen- 
dence of transport coefficients roughly follows the Arrhe- 
nius law. 

Schematic descriptions of the PES have often been in- 
voked to explain the fragile versus strong behavior of 
supercooled liquids ?i2, Strong glass- formers are expected 
to have a rough energy landscape, with energy barriers 
whose amplitude is essentially independent of the energy 
level. On the other hand, fragile glass-formers should 
display a more complex organization of stationary points 
and a broader distribution of energy barriers. Under- 
standing, at a quantitative level, the varying degree of 
fragility in different glass-formers represents a formidable 
task for theories. Correlations between fragility and sta- 
tistical or vibrational properties of local minima of the 
PES^ii^ have recently received a critical assessment for 
a wide range of models of the PES.^^ Variations in the 
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properties of the PES explored by supercooled liquids at 
different densities have also been discussed,— but their 
correlation to fragility have been a matter of debateJ^ 
Detailed studies have focused on the role of elementary 
rearrangements between adjacent local minima through 
transition states, both at constant densitj*i^ and con- 
stant pressurCi^^ Here we will concentrate our attention 
on the properties of high-order stationary points of the 
PES, whose relevance for supercooled liquids has been 
emphasized in the last years j-'^'^d^d^i^Oi^^i^^i^'^ High-order 
stationary points could offer a simple explanation of the 
fragile behavior of glass- formers, in terms of an increase 
of average energy barriersi^ This feature is encoded, 
in an effective way, in a number of models of energy 
landscapes developed in the last years r^i^^i^ZiiS^ and has 
sometimes been addressed in numerical simulations!^ 

Statistical properties of high-order stationary points of 
the PES have been investigated recently for a variety of 
monoatomic and binary systems, both in the liquid^^i^ 
and supercooled regime.— *2S The existence of some uni- 
versal features in the energy landscape of different model 
liquidsiii^ has been highlighted. At least in the case 
of the modified soft-sphere mixtures studied in Reffs^. 
such a universality has been found to refiect a fragility 
invariance of the systems investigated.-^S. In this work, 
we take a complementary point of view and ask: Are 
there variations in the properties of high-order station- 
ary points which correlate to fragility? To address this 
point, we consider a set of Lennard-Jones mixtures cooled 
at constant pressure: a series of equimolar, additive mix- 
tures with varying size ratio, together with some well- 
studied binary mixtures (Sec. |ll|. By investigating the 
temperature dependence of effective activation energies 
for relaxation (Sec. lIIl| . we provide support to our previ- 
ous results,-3^ which indicated the existence of systematic 
variations of isobaric fragility in additive mixtures and a 
remarkable pressure invariance in the mixture of Kob and 
Anderseni^i^ These trends allow us to test the connec- 
tion between fragility and some statistical properties of 
stationary points of the PES (Sec. lIVp . In particular, we 
show how fragility can be reflected in the saddles' den- 
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TABLE L Summary of thermal histories and simulation de- 
tails. Also shown are the number concentration of large par- 
ticles xi, the cutoff scheme used (see text for definitions) 
and the value of the cutoff radius r^- In the case of AMLJ- 
A mixtures, the following values of A have been considered: 
A = 0.60, 0.64, 0.70, 0.73, 0.76, 0.82. 
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sity of states, average energy barriers and localization 
properties of unstable modes. 



II. MODELS AND SIMULATION TECHNIQUES 

The binary mixtures studied in this work consist of 
N classical particles interacting via the Lennard- Jones 
potential 



(1) 



where a,/? — 1,2 are indexes of species. Each system is 
enclosed in a cubic box with periodic boundary condi- 
tions. In the following, reduced Lennard-Jones units will 
be used, i.e., cth, en and -y/ miali/eii as units of dis- 
tance, energy and time, respectively. Most simulations 
have been performed for samples of = 500 particles, 
employing the cutoff scheme of Stoddard and Ford^i at 
a cutoff radius Tc = 2.5. This cutoff scheme (QS) adds a 
shift and a quadratic term to the potential in order to en- 
sure continuity up to the first derivative of aX r — re- 
The role of the continuity of derivatives at the cutoff ra- 
dius Tc has been discussed in connection to minimization 
procedures 1^ To further investigate this point, we have 
also tested cut-and-shifted (CS) and cubic-splined cutoff 
(CSPL)^ on smaller samples composed by iV = 108 par- 
ticles. In this slightly smaller value of the cutoff 
radius has been used [tc = 2.2). 



1.8 T- 
1.7^ 
1.6, 
1.5^ 
1.4^ 
=^1.3- 
1.2- 
1.1 '- 
1.0^ 
0.9^ 
0.8 



AMLJ-0.64 
AMLJ-0.70 
AMLJ-0.76 
AMLJ-0.82 



0.0 0.5 1.0 1.5 



2.0 2.5 3.0 3.5 4.0 4.5 
T 



V 




□ P=5 
o P=10 
A P=20 









FIG. 1: Temperature dependence of density p(T) along iso- 
baric quenches at P = 10 for a selection of AMLJ-A mixtures. 
From bottom to top: A = 0.82, 0.76, 0.70, 0.64. 
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FIG. 2: Temperature dependence of density p{T) along iso- 
baric quenches at for BMLJ at P = 5, P = 10 and P = 20. 



We will focus our attention on the following binary 
Lennard-Jones mixtures: (i) The BMLJ mixture of Kob 
and Andersen, probably the most widely employed 
model for numerical simulations of the glass-transition, 
(ii) The WAHN mixture of Wahnstrom^^S which is an- 
other well-studied model glass-former, (iii) A set of ad- 
ditive, equimolar mixtures called AMLJ-A, characterized 
by different values of size ratio A = cr22/o'ii. In this case, 
the size ratio is varied in the range 0.60 < A < 0.82. A 
summary of all interaction parameters, together with a 
more detailed description of these models, can be found 
in Ref S 

Molecular Dynamics simulations were performed by 
cooling the liquid at constant pressure using Berendsen 
thermostat and barostat during the equilibration phase. 
Standard Velocity- Verlet algorithm was used to integrate 
the equations of motion. In order to achieve better con- 
trol on temperature in the deeply supercooled regime. 



3 



we performed a few production runs using the Nose- 
Poincare thermostat The timestep St was varied be- 
tween 0.002 at high temperature and 0.006 at low tem- 
perature. The time constant for the Berendsen thermo- 
stat^ was tT = St/ 0.1, while the coupling constant the 
for Berendsen barostat^ was 10'^ in reduced units. The 
inertia parameter of the Nose-Poincare thermostat was 
set to Q = 5. Constant pressure simulations provide a 
means to compare different mixtures in a way similar to 
the one employed in experiments. There is also interest in 
understanding how the sampling of the energy landscape 
changes when isobaric quenches are considered instead 
of isochoric ones.- ^ The density variations along isobaric 
quenches at a pressure P = 10 are shown for a selection 
of mixtures in Figs. [T] and [21 In order to make a com- 
parison with standard constant density simulations, we 
also performed some isochoric quenches for BMLJ and 
WAHN by fixing the density at the value used in the 
original papers. A summary of our thermal histories is 
shown in Table [J_For further details about quenching 
protocols see Refl34l. 

Description of the minimization technique employed 
for locating stationary points of the PES requires some 
further comments. We followed a simple and popular 
approachfii which consists of minimization of the mean 
square total force W of the system, using the L-BFGS 
algorithm,—. For each state point, some hundreds of inde- 
pendent configurations (typically between 200 and 1000) 
from simulation runs were considered as starting points 
for M^- minimizations. Some care has to be taken, since 
this numerical procedure often leads to quasisaddles, i.e., 
points with small but non-zero W, which display an 
inflection modcj ^^'^^i^^ As a criterion for distinguishing 
true saddles we use W < Wq = 10~^^, similarly to previ- 
ous studies jiii^ The fraction of true saddles sampled in 
small-sized samples {N — 108) is rather large (from 10% 
to 30% for WAHN, from 5% to 20% for BMLJ, depend- 
ing on temperature) , so that this approach appears to be 
quite feasible for similar system sizes. For larger samples, 
the fraction of true saddles decreases. However, from 
on overall point of view, our findings indicate that sad- 
dles and quasisaddles share similar average properties, in 
agreement with several previous observations ] ^^'^^i"^^ We 
also found that ensuring continuity of the interaction po- 
tential up to the first derivative at r = rc, i.e., employing 
the QS cutoff, is enough to avoid major round-off errors 
in M^-minimizations. Moreover, the fraction of true sad- 
dles sampled for TV = 108 does not increase appreciably 
when using the smoother CSPL cutoff. 



III. EFFECTIVE ACTIVATION ENERGIES 

A common way to display the temperature dependence 
of structural relaxation times t(T) in supercooled liq- 
uids is by construction of the so-called Angell plot, in 
which logr is shown against Tg/T. In a previous workf21 
we used a similar approach to analyze the variations of 



fragility in Lennard- Jones mixtures. Here, we take a dif- 
ferent view of the same problem and analyze the effective 
activation (free) energy E{T) defined by inversion of 



r(T) = Too exp 



EiT) 



T 



(2) 



where Tqo is the high-temperature limit of relaxation 
times. Analysis of the temperature dependence of E{T) 
will allow us to make contact with previous work based on 
effective activation energies;^ and to discuss the varia- 
tion of fragility in a way more convenient for highlighting 
the role of the PES (Sec. |IV|. 

Experimental and numerical analysis of effective ac- 
tivation energies E(T) in supercooled liquids clearly 
signals a crossover between two distinct dynamical 
regimes. At high temperature, relaxation times fol- 
low a mild, Arrhenius-type temperature dependence. 
Too exp [Eoo/T]. Hence, in the normal liquid regime, we 
have E{T) « Eoo- Below some crossover temperature 
T*, effective activation energies of fragile glass- formers 
increase markedly, indicating super- Arrhenius behavior. 
The more fragile is the glass-former, the shaper the in- 
crease of effective activation energies below T*. The 
features of E{T) just discussed are incorporated in the 
so-called frustration-limited domains theory of the glass- 
transitioufii which predicts E{T) to be of the form 



E{T) 



BT*(l- 



,8/3 



(3) 



It has been shown^ that the functional form in Eq. ([3]) 
provides a fair account of a wide spectrum of experimen- 
tal data. Fragility is measured by B, which is the param- 
eter quantifying the departure from the high-temperature 
Arrhenius regime. Further discussions on the role of the 
other parameters in Eg. (|3|) and on the exponent 8 /3 can 
be found in RefsH and47l. 

The dynamical quantity on which we focus in this sec- 
tion is the relaxation time for the decay of density fluc- 
tuations, as identified by the self part of the intermediate 
scattering function F"(k, t), where a — 1, 2 is an index of 
species. We define relaxation times Tq by requiring that 
F"{k*,t) has decayed to 1/e, where fc* w 8 is the posi- 
tion of the first maximum in the number-number static 
structure factor of the mixtures in consideration!^ A first 
guess of the crossover temperature T* is provided by the 
temperature Tonset at which two-step, non-exponential 
relaxation of F"{k,t) is first observed upon cooling the 
liquid from high temperature.=3S- For fitting our data to 
Eq. ([3]), we proceed as suggested by Kivelson et al.^ 
First we fit the high temperature portion of our data 
(T > Tonset) to an Arrhenius law Tqo exp [i?oo/T] and 
then we use Tqo and i?oo as fixed parameters in global a 
fit of our simulation data to Eq. ([3]). Note that, although 
T* is considered as a fitting parameter, its final value is 
never too far from the initial guess Tonset- 

In Table [Til we summarize the results of our fitting pro- 
cedure for Eq. ([3]). Considering separately the cases of 
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TABLE II: Parameters of fits to Eq. ([Sjl for effective activation energies E{T) of large and small species. The reference 
temperature T,. and the onset temperature Tonset are described in the text. 
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effective activation energies for large and small particles, 
we find that the fitted parameters for the two species 
show similar trends of variation in different systems and 
for different pressure and density conditions. In the fol- 
lowing, we will thus simply focus on the effective acti- 
vation energies E(T) obtained from the relaxation times 
T = Ti of large particles. In order to put into evidence 
the variation of fragility index B for different mixtures we 
plot, as in RefUJ, the difference {E{T)-E^)/T* against 
the reduced temperature T/T* . In Fig. [31 we show results 
obtained along isobaric quenches at P = 10 for a selec- 
tion of AMLJ-A mixtures (upper plot), and for BMLJ and 
WAHN mixtures (lower plot). The first important point 
is that there is a systematic variation of fragility with 
size ratio. Below the crossover temperature T* effective 
activation energies increase faster as the size ratio A in- 
creases, i.e., AMLJ-A mixtures become more fragile as A 
increases. The second point is that the BMLJ mixture is 
less fragile than the WAHN mixture. These observations 
are substantiated by the outcome of our fitting proce- 
dure. From an overall point of view, we find that Eq. ^ 
provides a good fitting function for our simulation data. 
Actually, the crossover around T* in our simulation data 
is smoother than predicted by Eq. ([3]), but it should also 
be remarked that Eq. ([S]) is not expected to hold exactly 
around T*4S, In the inset of the upper plot of Fig. [31 the 
fragility parameter B is shown as a function of size ra- 
tio for AMLJ-A mixtures. Despite the somewhat large 
uncertainty on our estimate of B, there is a clear trend 
of increase of B as A increases and a tendency to satu- 
rate around A « 0.80. Results obtained along different 
isobars for BMLJ show that the isobaric fragility index 
B for this system is essentially pressure invariant in the 
range 5 < P < 50, as it can be seen from the inset of the 



lower plot in Fig. [31 

Within the frustration-limited domains theory the 
crossover temperature T* is interpreted as an intrin- 
sic ridge between two distinct dynamical regimes, and 
should thus provide a means to scale and compare prop- 
erties of different glass-formers. Nonetheless, the use 
of isochronic, conventional reference temperatures, such 
as the glass-transition temperature Tg, is often useful 
and effective!^ We thus introduce a reference temper- 
ature Tr such that the relaxation time for large particles 
reaches T{Tr) = 4 x 10^. Since the value Tqo obtained 
from the high-temperature behavior is roughly system 
independent at a given pressure, the activation ener- 
gies for different systems converge to a common value 
E{Tr)/Tr « log(T(rr)/Too) « 12 as T ^ Tr. Thus, a 
plot in which both E{T) and T are scaled by Tr can 
be regarded as a generalized Angell plot, in which ac- 
tivation energies for different systems, when considered 
along the same isobar, converge around Tr to a common 
value. This kind of plot is shown in Fig. [H where we 
compare AMLJ-A mixtures for different values of A (up- 
per plot) and the two prototypical mixtures BMLJ and 
WAHN (lower plot). In this representation, fragility can 
be seen from a more pronounced increase of effective acti- 
vation energies, relative to the high temperature limit. A 
rough estimate of fragility can be thus obtained from the 
value of Eao/Tr, a fact which resembles the experimental 
correlation between Eoo/Tg and fragility 4^ 

A comparative analysis, based on Eq. ([3]), of exper- 
imental and numerical data was attempted some years 
ago by Ferrer et al.^ and further discussed by Tarjus et 
al.^ The outcome of the fitting procedure led these au- 
thors to raise some doubts about the fragile nature of 
some simulated models of supercooled liquids, including 
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FIG. 3: Effective activation energies for relaxation of large 
particles, after removal of the high-temperature limit Eoo- 
The dependence of {E{T) — Eoo)/T* on reduced temperature 
T /T* is shown along isobaric quenches at P = 10. Dashed 
lines are fits to Eq. Upper plot: AMLJ-A mixtures for 
values of size ratio A = 0.60 (squares) ,0.70 (triangles) and 
0.76 (circles). Lower plot: BMLJ (filled circles) and WAHN 
(open symbols). Inset of upper plot: fragility index B of 
AMLJ-A versus A. Inset of lower plot: fragility index B versus 
P obtained for different isobars in BMLJ. 



FIG. 4: Scaled effective activation energies E/Tr as a func- 
tion of T/Tr, along the isobar P — 10. Dashed lines are 
fits to Eq. dSl). Upper plot: AMLJ-A mixtures, for A = 0.60 
(squares), A = 0.70 (triangles) and A = 0.76 (circles). Lower 
plot: BMLJ (filled circles) and WAHN (open circles). 



the BMLJ mixture. This was contrary to the expecta- 
tion, based on qualitative grounds;^ that Lennard-Jones 
mixtures should be fragile glass-formers. Given the va- 
riety of Lennard-Jones models and external conditions 
analyzed in this work, we are probably in the position to 
shed some light on this point. First, we note that, for all 
mixtures considered, the ratio E{T)/Eoo is already larger 
than 2 around T^. We remark that this is a typical frag- 
ile behavior, even when compared to experimental data 
for fragile glass-formers such as ortho-terphenyl4i Note 
that, a part from the trivial determination of Eoo, this 
results is independent on the fitting procedure. Second, 
comparisons between experiments and numerical simula- 
tions of supercooled liquids should always be made with 
care. A much more limited temperature range is available 
in numerical simulations, and this can bias the results of 
fits to Eq. ([3]). For instance, by restricting the temper- 



ature range for fitting so that r < 10^, we obtained for 
BMLJ values of fragility index as low as i? « 12 at con- 
stant pressure, and i? w 4 at constant density, in line 
with the results obtained in RefH^ by considering a sim- 
ilar range of t. Fitting our data down to Tj., we obtain 
£? w 30 for BMLJ at constant pressure, and we expect 
that equilibrating the system at even lower temperatures 
would yield slightly larger values of B. Also note that 
for additive Lennard-Jones mixtures with moderate size 
asymmetry we find B w 100, which is already typical 
of intermediately fragile liquids {B w 90 for glycerol) 41 
Thus, from an overall point of view, Lennard-Jones mix- 
tures appear to be fragile glass-formers, as may be ex- 
pected for simple systems with non-directional interac- 
tions. On the other hand, it is true that some Lennard- 
Jones mixtures are less fragile than others. In particular, 
the well-studied BMLJ mixture, is not among the most 
fragile Lennard-Jones mixtures. 
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IV. POTENTIAL ENERGY SURFACE 

An important role for understanding the dynamical 
features of supercooled liquids is played by the station- 
ary points of the PES and, more generally, by the neg- 
atively curved regions of the PES.— i^SiSiiS^ In this sec- 
tion, we will adopt a simple, non-topographic point of 
view, ignoring the connectivity of stationary points. Ap- 
proaches based on pathways connecting adjacent min- 
ima through transition states,— or transitions between 
metabasins®'^'^ have been recently developed and applied 
to some model supercooled liquids, but they require ex- 
pensive and complex numerical procedures. We will thus 
focus on some simple statistical features of the PES, and 
discuss their correlations to fragility in Lennard-Jones 
mixtures. 

In the following, we will investigate the local curvature 
of the PES by looking at the Hessian matrix Ti of the 
potential energy. Standard diagonalization of Ti. yields 
a set of 2iN modes with eigenvalues and eigenvec- 
tors e", where a — 1, . . . , 3iV is an index of mode and 
j = 1, ... ,7V is an index of particle. Modes are then 
classified as stable if is positive (real frequency), or 
unstable if is negative (imaginary frequency). For 
liquids, most of the relevant information is encoded in 
the unstable modes of the PES, whose analysis usually 
comes in two varieties^^. The first approach is referred 
to as Instantaneous Normal Modes (INM) analysis, and 
focuses on instantaneous configurations sampled along 
the MD trajectory i^Si^i^ The second approach consid- 
ers high-order stationary points of the PES, obtained 
using minimization procedures fii'^^'i^i'^'^ According to 
the number of unstable modes in the Hessian matrix, 
stationary points are classified as local minima {uu — 0) 
or saddles (n„ 7^ 0). As mentioned in Sec|lTl quasisad- 
dles are other points of the PES often reached by the 
minimization algorithm employed in this work. The ex- 
clusion of quasisaddles from statistical averages will not 
affect the main conclusions of this section. We checked 
the reliability of our results on some smaller samples of 
iV = 108 particles, in which a larger fraction of saddles 
could be found (Sec. [11} • In the following, we will focus 
on the larger sample [N = 500) and we will mostly use 
the term saddles in a broad sense, without distinction 
between saddles and quasisaddles. 

As a starting point, we consider the ensemble- averaged 
density of states (DOS) 



37V 



p(c.;T) = (^<5( 



- w, 



(4) 



at temperature T. The thermal average in Eq. ^ can be 
performed using either instantaneous configurations (i- 
DOS) or saddles (s-DOS). The unstable branch of p(w; T) 
will be denoted by pu{uJ',T), and imaginary frequencies 
will be shown, as usual, along the real negative axis. To 
provide a quantitative account of the features of Pu{^, T), 
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FIG. 5: Unstable branch of density of states Pu{u>) for in- 
stantaneous configurations (i-DOS, upper plot) and saddles 
(s-DOS, lower plot) in BMLJ. Results are shown at four dif- 
ferent state points at P = 10 (from top to bottom, T = 
2.0, 1.0, 0.7, 0.6). Dashed lines are fits to Eq. ©. 



we consider the following functional form 



Pui^J] T) ^ Alo exp 



(5) 



which has been shown to describe well the unstable i- 
DOS.^?. Specific functional forms for the s-DOS have been 
discussed, for instance, in the context of p-spin mod- 
els f"-''— but they fail to meet some basic requirements 
for realistic systems, such as the behavior p{uj) for 

— > 0. We have thus attempted to apply Eq. ^ also to 
s-DOS and found that Eq. ((5]) provides an excellent fit 
for both i-DOS and s-DOS, becoming inadequate only at 
very high temperature or in the limit of vanishing interval 
of imaginary frequencies. The quality of the fits obtained 
using Eq. ([5]) is exemplified in Fig. [5] for different state 
points of BMLJ. 

In the context of the INM theory of diffusion,— super- 
Arrhenius behavior is explained in terms of the temper- 
ature dependence of the parameters in Eq. ([5]) for the 
i-DOS, and should be primarily signaled by an increase 
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0.64 (filled circles) and AMLJ-0.82 (open circles) at P = 10. 
Lower plot: BMLJ at different pressures, P = 5 (squares), 
P = 10 (circles) and P = 20 (triangles). Dashed lines are 
fits of the type a + b{Tr/T), with 6 ~ 5 being roughly system 
independent. . 



of C by decreasing temperature i^^'^^i^^ '^ Unfortunately, 
we found that INM theory is not able to put into evidence 
the different fragility of the Lennard- Jones mixtures con- 
sidered in this workS. For instance, we found that C, 
as a function of T/Tr, has similar values in all systems, 
whereas we would have expected a sharper increase in 
the case of more fragile mixtures. Analysis of the s-DOS 
in terms of Eq. ^ provides a different, sharper picture. 



Without attempting a detailed analysis of the tempera- 
ture dependence of all parameters in Eq. ([S]), we will fo- 
cus on the behavior of parameter C. In Fig. [SI we show 
the dependence of C on T/Tr for the s-DOS of different 
mixtures at constant pressure. Apart from some devia- 
tions at very high-temperature, C decreases by decreas- 
ing temperature, differently from the case of the i-DOS. 
In all systems, we observe a temperature dependence of 
the type C ~ b{Tr/T). Interestingly, all mixtures seem 
to share a common value of the slope 6 « 5 in a plot of 
C versus Tr/T, and we find a shift towards larger val- 
ues of C as fragility increases. As a check of the relation 
between C and fragility, the values of C along different 
isobars for BMLJ collapse on a master curve when plot- 
ted against T/Tr, consistently with the pressure invari- 
ance of fragility reported in Sec. IIIII Our results thus 
indicate that Eq. ([5]) could provide a good starting point 
for modeling the s-DOS in realistic systems and that a 
saddle-based approach is more sensitive to the dynami- 
cal behavior of supercooled Lennard- Jones mixtures than 
INM. 

In the light of previous studies of vibrational properties 
of local minimafiiii^ it might be asked whether coarse- 
grained quantities obtained from the s-DOS, such as the 
average frequency of of stable modes 



uJsiT) 



and unstable modes 



dij Ljpiuj; T) 



duj ujp{uj; T) 



(6) 



(7) 



already convey information about fragility. In Fig. [71 we 
show the dependence of cou and on T/Tr for different 
mixtures at constant pressure. We found analogous ther- 
mal behaviors by considering geometric mean frequen- 
cies of stable and unstable modes. Similarly to what 
happens in the case of local minima, constant-pressure 
data show an increase of lOs by decreasing temperature, 
i.e., by decreasing energy of saddles. This behavior is 
opposite to the one observed in constant-density simu- 
lations. The average frequency of unstable modes uJu 
always shows a non-monotonic temperature dependence, 
characterized by a maximum at intermediate tempera- 
tures, which is peculiar to isobaric quenches. Comparing 
mixtures along isobaric quenches at P = 10, we find a 
slight shift to larger absolute vibrational frequencies, as 
fragility decreases. However, the robustness of this cor- 
relation is weakened when it is tested using the pressure 
invariance of isobaric fragility in BMLJ. In the bottom 
plot of Fig. [71 we look at the behavior of cug and uJu along 
different isobars in BMLJ. As the pressure P of the isobar 
increases, vibrational frequencies are shifted markedly to 
larger absolute values, most probably by the increasing 
density®^. This behavior led us to reconsider the case 
of local minima along different isobars in BMLJ, and we 
found a similar trend in vibrational properties. Thus, 
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FIG. 7: Average frequency of stable modes los (main plots) 
and unstable modes uJu (insets) of saddles as a function of 
T/Tr. Upper plot: BMLJ (filled circles) and WAHN (open 
circles) at P = 10. Middle plot: AMLJ-0.64 (filled circles) 
and AMLJ-0.82 (open circles) at P = 10. Lower plot: BMLJ 
at P = 5 (squares), P = 10 (circles), and P = 20 (triangles). 



although some correlation might be observed at a given 
pressure, there seems to be no direct connection between 
average vibrational frequencies of stationary points and 
fragility. 

The variation of fragility in Lennard- Jones mixtures, 
as discussed in terms of effective activation energies for 
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FIG. 8: Upper plot: scatter plot of the fraction of unstable 
modes against energy of single saddles. Results are shown 
for WAHN at P = 10 for three different state points: T = 
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estimate the derivative in Eq. (Ilip . i.e., Es — l/3b. Lower 
plot: parametric plot of average unstable modes of saddles 
fu{T) against energy of saddles es(T), for WAHN at P = 10. 



relaxation E{T), calls for an explanation based on energy 
barriers. Whereas it is clear that E{T) in Eq. ^ is 
rather an activation free energy, the leading contribution 
to it might already come from potential energy barriers 
between single saddles. To address this point, we follow 
the simple proposal of Cavagnai^ The starting point is 
the relation 



(8) 



between the fraction of unstable modes and the energy 
of saddles. Eq. ([5]) will be treated as parametric in T, 
i.e., we consider the average fraction of unstable modes 

fu^fu{T) = {nu/m)^= [ dujpu{cu;T) (9) 



and the average energy of saddles 

es = es(T) 



(10) 



at temperature T. It has been shown that Eq. ([8]), as 
obtained from numerical simulations, is insensitive to the 
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actual minimization algorithm employed,-^- and to the 
inclusion of quasisaddles.^? According to Cavagna~ the 
average energy difference 

Esies) = (11) 
3 df u 

between saddles of order n and n + 1 provides an esti- 
mate of potential energy barriers in the PES. More re- 
fined treatments would take into account the connectivity 
of saddles and existence of a distribution of energy barri- 
ers.^'* In order to find Es(T), we compute the derivative 
in Eq. pT|) by linear regression of Cs vs. scatter data 
of saddles sampled at temperature T, as illustrated in 
Fig.li 

The temperature dependence of the effective energy 
barriers Eg (T) is shown in the left plots of Fig. [9] for 
different mixtures at constant pressure. Below T*, i.e., 
in the range of temperature where activated dynamics is 
expected to become important fii^ the behavior of Eg (T) 
correlates to the fragility of the mixture. In fact, the in- 
crease of effective energy barriers upon supercooling is 
sharper and more pronounced, the more fragile is mix- 
ture. In the case of the more fragile mixtures, we find 



a striking similarity between the increase of Es{T) be- 
low T* and that of the effective activation energies E{T) 
defined by Eq. In WAHN, for instance, we find 

E{Tr) « Es{Tr) « UTr. The trends just discussed are 
in line with the results obtained by direct calculations 
of energy barriers between adjacent minima in the soft- 
sphere version of WAHNS^ and BMLJii Some concerns 
might regard the fact that es(T), i.e., the energies of 
saddles sampled at a given T, can depend on the mini- 
mization algorithm4i On the other hand, the results ob- 
tained in Ref.lii indicate that the energy dependence of 
the properties of saddles is much less sensitive to the de- 
tails of the minimization procedure employed. We have 
thus analyzed our data for Eg as a function of , where 
Cs is given by the thermal average in Eq. (jlOp . focusing 
on the energy range below es(T*). For convenience, we 
have shifted the energies Cs by es{T*). Such a represen- 
tation of our data is shown in the right plots of Fig. [§] 
and confirms the trends discussed above on the basis of 
the temperature dependence of Es. Thus, independent 
of the representation used, the average energy barriers 
show a strong connection to the variations of fragility in 
our models. This also provides evidence of the relevance, 
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for the supercooled dynamics, of the es(T) mapping ob- 
tained through VF-minimizations. 

What is the effect of pressure on energy barriers? From 
the plot in Fig. 1101 we see that increasing pressure in 
BMLJ leads to larger potential energy barriers. This 
behavior is consistent with the results obtained by Mid- 
dleton and Wales who calculated the distribution of 
potential energy barriers for diffusive rearrangements at 
different pressures for BMLJ. What is made clear by our 
results, is that, at least in the case of BMLJ, the in- 
crease of potential energy barriers with pressure has lit- 
tle dynamical impact, because it is compensated by the 
increase of the reference temperature T^. That is, larger 
energy barriers will be sampled at higher temperatures. 
Starting from data along different isobars in BMLJ, in 
fact, we could obtain a rough master curve by scaling 
both Es and T by T^. 

The unstable modes of saddles sampled in the super- 
cooled regime may provide information about the ele- 
mentary dynamical processes in the systemj^ It is thus 
of interest to analyze the spatial localization features of 
unstable modes in different Lennard- Jones mixtures, and 
see how they relate to fragility. To address this point, we 
first average the squared displacements for each particle 
over all the unstable modes^"^ 



ri u ^ 

Q — 1 



(12) 



Two different measures of localization for the vector of 
average squared displacements Ef^ are then considered. 
The reduced gyration radius is defined as 



1 



L/2 



■ N 

El 

.2=1 



1/2 



(13) 



where L is the side of the simulation box. This quan- 
tity equals 1 when the vector Ef"^ is extended over the 
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FIG. 11: Participation ratio of average squared displace- 
ments on unstable modes as a function of reduced tempera- 
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whole system, and decreases progressively as the spatial 
localization of Ef"^ becomes more pronounced. The par- 
ticipation ratio is defined as 



- N 



N 

E 

1=1 



Er 



(14) 



and provides a rough estimate of the fraction of particles 
having significant displacements in -E"^. For instance, 
should be 0(1) when the unstable modes are homo- 
geneously distributed in the system. The temperature 
dependence of these two quantities is shown in Fig. [TT] 
for different Lennard- Jones mixtures. The existence of a 
sharp localization of unstable modes around Tr, as iden- 
tified by the abrupt decrease of L'^{T), appears to be 
a universal feature of saddles sampled by supercooled 
Lennard- Jones mixtures. On the other hand, the pat- 
tern of localization of the unstable modes changes ac- 
cording to the fragility of the mixture. The more fragile 
is the mixture, the more rapid the localization of unstable 
modes upon supercooling, as it is suggested by the behav- 
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lor of p"(r). In the range of temperature above T^, we 
find that fragile mixtures tend to have larger values of p". 
In this case, a larger fraction of particles is thus involved 
in the unstable modes, which is consistent with expecta- 
tion that rearrangements should involve larger clusters as 
fragility increases?^ We found further support to these 
considerations by analyzing the average participation ra- 
tio and gyration radius of individual unstable modes of 
saddles. 

Inspection of animated unstable modes of saddles 
sampled at low temperature indicates the occurrence 
in BMLJ of strongly localized, high-frequency unstable 
modes, in which few small particles show very large dis- 
placements. This feature is reflected in a clearly bimodal 
distribution of i?"^ for small particles in deeply super- 
cooled BMLJ. In Fig. [T^l we show the distributions of 
^ for BMLJ and WAHN at the lowest equilibrated 



temperatures. In the case of BMLJ, in fact, we find a 
bump at large values in the distribution of E""^ for small 
particles. We also often observed correlated, string- like 
rearrangements of large particles in the unstable modes 
of BMLJ. This feature is exemplified in the snapshots 
of Fig. [T31 where we show the unstable modes of a typ- 
ical quasisaddle sampled in deeply supercooled BMLJ. 
By comparison, unstable modes in WAHN tend to in- 
volve larger and more compact clusters of particles and 
to possess a more pronounced spatial overlap. These 
features, in the light of our previous investigations^^ 
should influence the dynamical processes within the (3- 
relaxation timescale, and could provide the basis for un- 
derstanding the microscopic origin of dynamical hetero- 
geneities^ on longer timescales. Analysis of the connec- 
tivity between stationary points could also help explain- 
ing the relative weight of different dynamical processes — 
string-like motionsj^^'-- democratic rearrangements^^ — 
observed in supercooled Lennard- Jones mixtures. More 
detailed studies along this direction will require signifi- 
cant additional work. 



V. CONCLUSIONS 

Molecular Dynamics simulations of supercooled 
Lennard- Jones mixtures continue to provide a useful 
benchmark for theories and paradigms of the glass- 
transition. A minimal exploration of the field of param- 
eters of the Lennard-Jones potential for binary mixtures 
has revealed a rich phenomenology. In particular, we 
found a systematic variation of fragility, i.e., a varying 
degree of super-Arrhenius behavior of dynamical prop- 
erties. By analyzing the temperature dependence of ef- 
fective activation energies for relaxation, we found that 
fragility increases by increasing size ratio A — (722 /cn 
in equimolar, additive mixtures. Two prototypical mix- 
tures, the BMLJ mixture of Kob and Andersen^^^ and the 
one of Wahnstrom,— also have different fragility indexes. 
As an interesting result, we also found that fragility does 
not change appreciably with pressure in BMLJ. 

In a previous paper;^ we discussed these trends in 
terms of the thermal rate of growth of locally preferred 
structures upon supercooling. Here, we have investigated 
the different, complementary role of the potential energy 
surface explored in the supercooled regime. We have 
adopted a simple, non-topographic approach and ana- 
lyzed some statistical properties of the PES, with par- 
ticular focus on high-order stationary points. We have 
provided an estimate of average potential energy barri- 
ers and found a striking correlation with the fragility of 
the mixture: the more fragile the mixture, the more pro- 
nounced the increase of energy barriers upon supercool- 
ing. Not ignoring the role of metabasins^ and multistep 
processes,— an increase of energy barriers between sin- 
gle saddles appears to be a simple, possible mechanism 
for super-Arrhenius behavior of dynamical properties in 
fragile glass-formers. We have also found that a proper 
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FIG. 13: Selection of three unstable modes of a quasisaddle {riu = 4) sampled in BMLJ at T = 0.66, P = 10. The nearly-zero 
mode of the quasisaddle has been ignored. A fourth unstable mode, not shown, is very similar in extension and shape to that 
shown in (b). For clarity, only particles having square displacements (e^)^ larger than 0.004 are shown, and eigenvectors are 
scaled logarithmically. Large and small particles are shown as pale large spheres and small darker spheres, respectively. Note 
the strong localization of mode (a) and the existence of distinct string- like instabilities of large particles in mode (b) and (c). 



characterization of the saddles' density of states will al- 
ready encode the relevant information about fragility. On 
the other hand, mean frequencies of stable and unsta- 
ble modes do not provide robust correlations to fragility. 
Their strong variations with density along different iso- 
bars in BMLJ, in fact, are not accompanied by a signifi- 
cant change in fragility. 

As a general rule, unstable modes of saddles become 
more and more localized upon supercooling, but this fea- 
ture is sharper and more pronounced, the more fragile is 
the mixture. This can be interpreted as the counterpart 
of a more rapid growth, upon supercooling, of slow do- 
mains, characterized by distinct locally preferred struc- 
tures 4^ From some preliminary calculations, we have 
found, as expected, that particles at the center of locally 
preferred structures are stabilized and are not involved 
in the unstable modes. Thus, the study of the potential 
energy surface presented in this work and our previous 
microstructural analysis^ complement each other very 
well. Formation of stable, long-lived structures, such as 
icosahedral domains in additive mixtures, could corre- 
spond to deeper traps in the energy landscape, thus forc- 



ing relaxation over larger energy barriers. On the other 
hand, frustration of stable prismatic domains could be 
the origin of the limited growth of potential energy bar- 
riers in BMLJ. The two approaches together may thus 
provide an intriguing picture of the fragile vs. strong 
behavior of glass-former, bridging the ideas of frustra- 
tion in supercooled liquids^ and roughness of the energy 
landscaped Assessment of such picture by studying dif- 
ferent interactions remains an open problem for further 
investigations. 

Acknowledgments 

The authors would like to thank A. Cavagna for use- 
ful discussions and a critical reading of the manuscript. 
Computational resources for the present work have been 
partly obtained through a grant from "Iniziativa Trasver- 
sale di Calcolo Parallelo" of the Italian CNR-Istituto 
Nazionale per la Fisica della Materia (CNR-INFM) and 
partly within the agreement between the University of 
Trieste and the Consorzio Intenmiversitario CINECA 
(Italy). 



^ F. Sciortino, J. Stat. Mech.: Theory Exp. p. P05015 
(2005). 

^ D. J. Wales, Energy Landscapes (Cambridge University 

Press, Cambridge, 2003). 
^ F. H. Stillinger and T. A. Weber, Phys. Rev A 25, 983 

(1982). 

* S. Sastry, P. G. Debenedetti, and F. H. Stillinger, Nature 

393, 554 (1998). 
^ F. H. Stillinger, Phys. Rev. B 41, 2409 (1990). 
^ B. Doliwa and A. Heuer, Phys. Rev. Lett. 91, 235501 

(2003). 



^ B. Doliwa and A. Heuer, Phys. Rev. E 67, 031506 (2003). 
® B. Doliwa and A. Heuer, Phys. Rev. E 67, 030501 (2003). 
^ R. A. Denny, D. R. Reichman, and J.-P. Bouchaud, Phys. 

Rev. Lett. 90, 025503 (2003). 
^° F. H. Stillinger, Science 267, 1935 (1995). 
" S. Sastry, Nature 409, 164 (2001). 

R. J. Speedy, J. Phys. Chem. B 103, 4060 (1999). 

G. Ruocco, F. Sciortino, F. Zamponi, C. De Michele, and 

T. Scopigno, J. Chem. Phys. 120, 10666 (2004). 

G. Tarjus, D. Kivelson, S. Mossa, and C. Alba-Simionesco, 

J. Chem. Phys. 120, 6135 (2004). 



13 



T. F. Middleton and D. J. Wales, Phys. Rev. B 64, 024205 

(2001) . 

T. F. Middleton and D. J. Wales, J. Chem. Phys. 118, 
4583 (2003). 

L. Angelani, G. Ruocco, M. Sampoli, and F. Sciortino, J. 
Chem. Phys. 119, 2120 (2003). 
1* J. P. K. Doye and D. J. Wales, J. Chem. Phys. 116, 3777 

(2002) . 

^3 D. J. Wales and J. P. K. Doye, J. Chem. Phys. 119, 12409 

(2003) . 

M. Sampoli, P. Benassi, R. Eramo, L. Angelani, and 

G. Ruocco, J. Phys.: Condens. Matter 15, S1227 (2003). 

L. Angelani, R. Di Leonardo, G. Ruocco, A. Scala, and 

F. Sciortino, Phys. Rev. Lett. 8, 5356 (2000). 

L. Angelani, R. Di Leonardo, G. Ruocco, A. Scala, and 

F. Sciortino, J. Chem. Phys. 116, 10297 (2002). 

D. Coslovich and G. Pastore, Europhys. Lett. 75, 784 

(2006) . 

A. Cavagna, Europhys. Lett. 53, 490 (2001). 

F. Zamponi, L. Angelani, L. F. Cugliandolo, J. Kurchan, 
and G. Ruocco, J. Phys. A: Math. Gen. 36, 8565 (2003). 
T. Keyes, J. Chowdhary, and J. Kim, Phys. Rev. E 66, 
051110 (2002). 

A. Andronico, L. Angelani, G. Ruocco, and F. Zamponi, 
Phys. Rev. E 70, 041101 (2004). 

A. Cavagna, I. Giardina, and G. Parisi, J. Phys. A: Math. 
Gen. 34, 5317 (2001). 

T. Grigera, A. Cavagna, L Giardina, and G. Parisi, Phys. 
Rev. Lett. 88, 055502 (2002). 
^° P. Shah and C. Chakravarty, J. Chem. Physics 115, 8784 
(2001). 

S. N. Chakraborty and C. Chakravarty, J. Chem. Phys. 
124, 014507 (2006). 

L. Angelani, C. De Michele, G. Ruocco, and F. Sciortino, 
J. Chem. Phys. 121, 7533 (2004). 

C. De Michele, F. Sciortino, and A. Coniglio, J. Phys.: 
Condens. Matter 16, L489 (2004). 

D. Coslovich and G. Pastore, J. Chem. Phys. 127, 124504 

(2007) . 

W. Kob and H. C. Andersen, Phys. Rev. E 51, 4626 (1995). 
W. Kob and H. C. Andersen, Phys. Rev. E 52, 4134 (1995). 
S. D. Stoddard and J. Ford, Phys. Rev. A 8, 1504 (1973). 
3® P. Shah and C. Chakravarty, J. Chem. Phys. 118, 2342 
(2003). 

G. Wahnstrom, Phys. Rev. A 44, 3752 (1991). 

*° S. D. Bond, B. J. Leimkuhler, and B. B. Laird, J. Comp. 
Phys. 151, 114 (1999). 
S. Nose, J. Phys. Soc. Jap. 70, 75 (2001). 
M. P. Allen and D. J. Tildesley, Computer Simulation of 
Liquids (Clarendon Press, 1987). 



D. C. Liu and J. Nocedal, Math. Program. 45, 503 (1989). 
T. S. Grigera, J. Chem. Phys. 124, 064502 (2006). 
"•^ D. Kivelson, G. Tarjus, X. Zhao, and S. A. Kivelson, Phys. 
Rev. E 53, 751 (1996). 

M. L. Ferrer, H. Sakai, D. Kivelson, and C. Alba- 
Simionesco, J. Phys. Chem. B 103, 4191 (1999). 
G. Tarjus, D. Kivelson, and P. Viot, J. Phys.: Condens. 
Matter 12, 6497 (2000). 

V. N. Novikov, Y. Ding, and A. P. Sokolov, Phys. Rev. E 
71, 061501 (2005). 

C. A. Angell, B. E. Richards, and V. Velikov, J. Phys.: 
Condens. Matter 11, A75 (1999). 

S. D. Bembenek and B. B. Laird, J. Chem. Phys. 104, 
5199 (1996). 

J. D. Gezelter, E. Rabani, and B. J. Berne, J. Chem. Phys. 
107, 4618 (1997). 

T. Keyes, G. V. Vijayadamodar, and U. Zurcher, J. Chem. 
Phys. 106, 4651 (1997). 

K. Broderix, K. K. Bhattacharya, A. Cavagna, A. Zip- 
pelius, and L Giardina, Phys. Rev. Lett. 85, 5360 (2000). 
^'^ A. Cavagna, L Giardina, and T. Grigera, J. Phys. A: Math. 
Gen. 36, 10721 (2003). 

T. Keyes, J. Chem. Phys. 101, 5081 (1994). 
"'^ G. Vijayadamodar and A. Nitzan, J. Chem. Phys. 103, 
2169 (1995). 

W.-X. Li and T. Keyes, J. Chem. Phys. Ill, 5503 (1999). 
'^^ L. Berthier and J. P. Garrahan, Phys. Rev. E 68, 041201 
(2003). 

^'^ E. A. Jagla, Mol. Phys. 99, 753 (2001). 
®° A. Widmer-Cooper, P. Harrowell, and H. Fynewever, Phys. 
Rev. Lett. 93, 135701 (2004). 

C. Donati, S. C. Glotzer, P. H. Poole, S. J. Plimpton, and 
W. Kob, Phys. Rev. E 60, 3107 (1999). 
T. B. Schroeder, S. Sastry, J. C. Dyre, and S. C. Glotzer, 
J. Chem. Phys. 112, 9834 (2000). 

G. A. Appignanesi, J. A. Rodriguez Fris, R. A. Montani, 
and W. Kob, Phys. Rev. Lett. 96, 057801 (2006). 
G. Tarjus, S. A. Kivelson, Z. Nussinov, and P. Viot, J. 
Phys.: Condens. Matter 17, R1182 (2005). 

^'^ J. Chowdhary and T. Keyes, Physica A 314, 575 (2002). 
Actually, also a third way has been considered.— 
Introduction of a lower frequency cutoff ujc to filter some 
shoulder modes^^ would not affect our conclusions. 

®* It has been suggestectiS. that the leading contribution to 
the density dependence of the geometric mean frequency 
in local minima should scale as a power law of p. We argue 
that a similar argument might hold for high-order station- 
ary points. 



